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WENO schemes for multidimensional 
nonlinear degenerate parabolic PDEs 
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Abstract 


In this paper, a scheme is presented for approximating solutions of non- 
linear degenerate parabolic equations which may contain discontinuous solu- 
tions. In the one-dimensional case, following the idea of the local discontinu- 
ous Galerkin method, first the degenerate parabolic equation is considered as 
a nonlinear system of first order equations, and then this system is solved us- 
ing a fifth-order finite difference weighted essentially nonoscillatory (WENO) 
method for conservation laws. This is the first time that the minmod-limiter 
combined with weighted essentially nonoscillatory procedure has been applied 
to the degenerate parabolic equations. Also, it is necessary to mention that 
the new scheme has fifth-order accuracy in smooth regions and second-order 
accuracy near singularities. The accuracy, robustness, and high-resolution 
properties of the new scheme are demonstrated in a variety of multidimen- 
sional problems. 


Keywords: WENO schemes; Finite difference scheme; Multidimensional 
nonlinear degenerate parabolic equation; Porous medium equation. 


1 Introduction 


This paper is concerned with multidimensional of nonlinear degenerate 
parabolic equations of the form 


ut = Ab(u) + S(x,t,u), x =(%,...,%¢) € Nand t>0, (1) 


satisfying the initial condition 
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u(x, 0) aa uo(X), XE Q, (2) 


and suitable boundary conditions. The condition Dj(u) = bj(u) = 0, 
j =1,...,d is needed to make the equation formally parabolic. Whenever 
D,(u) = 0, for some u € R, is said that the equation degenerates at that 
u—level, since it ceases to be strictly parabolic. For b;(u) =u, m > 1, and 
S(x,t, u) = 0, the equation ({) is called porous medium equation(PME) [B6]. 
As is said in |B6], if the initial condition has compact support, then (f) and 
(2) have at most one weak solution that has compact support in space (for 
more detail, see Theorem 5.3 of [R6]). The degenerate parabolic equations 
appear often in applications, for example, spread of viscous fluids []], math- 
ematical biology [7], groundwater [6], heat transfer [16], flow of an isentropic 
gas through a porous medium [86], and other fields. 


Methods for approximating solutions to (MI) have attracted a lot of atten- 
tion in recent years, for example, local discontinuous Galerkin finite element 
method [B9], relaxation scheme [9], and kinetic scheme [4] have been pro- 
posed to find its solution. Also, Liu, Shu, and Zhang [26] proposed two 
different formulations of WENO schemes for solving nonlinear degenerate 
parabolic equations. The first formulation approximates directly the second 
derivative term by using a conservative flux difference, while the second for- 
mulation is similar to the formulation of the local discontinuous Galerkin 
(LDG) schemes [1]. 


As we know, WENO schemes were introduced in the literature to approxi- 
mate hyperbolic conservation laws. Also, WENO schemes are based upon the 
essentially nonoscillatory (ENO) schemes [IZ]. The first version of WENO 
schemes was introduced in finite volume formulation for the one-dimensional 
conservation laws in 1994 by Liu, Osher, and Chan [25]. In 1996, Jiang 
and Shu observed that the ENO stencil selection is sensitive to the round-off 
perturbation near zeros of the solution and its derivatives [ZI]. Then, they 
improved the ENO reconstruction and proposed a general framework for de- 
signing arbitrary order accurate finite difference WENO schemes, which are 
more efficient for multidimensional calculations. 


The degenerate parabolic equation (Ml) may have discontinuous solution, 
possible existence of sharp fronts, and finite speed of propagation of wave 
fronts [B6]. The degenerate parabolic equations usually have features similar 
to those of a hyperbolic conservation laws. Therefore for solving (), it is rea- 
sonable to generalize numerical methods for solving hyperbolic conservation 
laws, such as the WENO schemes [26]. 


In this paper, a WENO finite difference scheme is proposed by following 
the idea of the local discontinuous Galerkin method [10]. First rewrite (1) in 
one-dimensional case as a nonlinear system of first order equations 


Up = Vz + S(ax, tu) 
v = d(u)a, 
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and then to solve the two first order equations using a fifth order finite differ- 
ence WENO method for conservation laws. This approach has the advantage 
of simplicity, and it is easier to generalize it to higher than second order 
PDEs. 


The structure of this paper is as follows: In section 2, we describe by 
detailing the construction and implementation of the scheme, for nonlinear 
degenerate parabolic equations. In section 3, the results of numerical exper- 
iments conducted with the proposed scheme are given. Some remarks are 
made in section 4. 


2 Nonoscillatory reconstruction to the second derivative 


In this section, a fifth-order WENO finite difference scheme for multidimen- 
sional nonlinear degenerate parabolic equation is described. Using the di- 
mension by dimension approach [2], the presented scheme in this paper for 
multidimensional problems is implemented. 


2.1 General framework 


Consider the one-dimensional nonlinear degenerate parabolic equation 


uz = b(U)oe + S(a,t, u), (a, t) € Q x (0,00), (3) 
u(a,0) = uo(z), 


satisfying suitable boundary conditions. 


This section describe the formulation of a WENO finite difference scheme 
by following the idea of the local discontinuous Galerkin method [1]. First, 
the equation (B) is considered as a nonlinear system of first order equations 


Ut = Uz + S(z, t, u) (4) 
v= 0(u)s, (5) 


and then to solve the two first order equations, respectively, using the fifth or- 
der finite difference WENO method for conservation laws [I]. This approach 
has the advantage of simplicity, and it is easy to generalize it to higher than 
second order PDEs. The effective stencil, which is a composition of two suc- 
cessive WENO procedures, is also wider in comparison with the approach 
that is said in [2]. To make the final effective stencil smaller, a right-biased 
stencil for (@) followed by a left-biased stencil for (B) is used. 
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2.2 The fifth order finite difference WENO method 


Assuming that the cell averages u; are known at time t, van Leer [85] proposed 
a piecewise polynomial reconstruction u(x,t) = )), Rj(x)x;(x) to solve 1- 
dimension hyperbolic conservation laws, where xj;(x) is the characteristic 


function of the cell Ij = [w;_1,2;41]. Let U(x) be the primitive function of 


u(a,t); that is, U(2) = ee u(€,t)d&. To begin, an optimal polynomial is 
selected, which is denoted by popt,;(x), to approximate u(x,t) on cell I; of 
degree four on the five-cells central stencil S = {I;-2,...,Ij;+2}. In [BO,BQ), 
this optimum polynomial is obtained. First, U(a) is interpolated at the 
cell boundaries {js penny Dips } by using Newton’s interpolation formula. 
Afterwards, differentiating the new polynomial, we get 


5 


t—-1 i—1 
Popt,j(@) = So Ula;-s, _ T5448] ‘> II (x — Byaj8)) 


i=l m=01=0,l4m 


where U[---] is a divided difference of the function U(x). High-order ENO 
reconstructions decrease damping of the solution, but generate significant 
oscillations when solving hyperbolic systems unless costly characteristic de- 
compositions are used. For solving this problem, three supplementary poly- 
nomials are defined, {p)(x), p{(x), and p3(x)}, approximating u(x) with a 
lower accuracy on J;. The three-cells stencils of ENO3 are obtained by either 
choosing r cells to the left and s cells to the right of J;. Then, the three-cells 
stencils of ENO3 can be explicitly determined. 


S” = {Ij-r,...,Tj4s}, (r,s) = (0,2), (1,1), and (2,0). 


Having obtained pop(x) as above, p'")4 (x) is obtained from cell boundaries 
of the stencil S” 


3 i-1 i-l 
O'G)=) Uta swwtwal>, Tl @-jaaa). © 
$= m=01=0,l4Am 


In order to obtain a more efficient reconstruction, the ENO3 reconstruction is 
modified. We interpolate equation (@) over an additional point lying within 
the same stencil and the same number of selection steps as ENO3. The new 
reconstruction starts by seeking a new high-order interpolating polynomial, 


p.(a), such that equation (@) passes over the additional point x; € J;, 
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The divided difference U[x;_,+8,25], for r = 0,1, and 2 of (@), is given by 
Ue; gst) = soe [JT wleas — Jc"? wea] 
j-rt+3 


er 
= re oe mg 63, L;(2)x i(2)) der. 
As it can be seen from the equation (8), a polynomial that retains infor- 
mation within the cell I; is needed. The NT scheme [29] is based on the 
first-order Lax—Friedrichs scheme [13] and it involves the reconstruction of 
piecewise-linear MUSCL-type interpolates from piecewise constant data and 
uses nonlinear limiters to prevent oscillations 


1 
3) Uj, ze I;. (9) 
Integrating over (Q), the divided differences in (Q) are given by 


1 
5 — 2r 


U|aj-r48 05] = («; (r= 2)tat 


1 
(r? _ 3r + 2)Uj+2 + i): r= 0, 1, 2. 


For the numerical derivative wu’, the UNO limiter [[§] is chosen which is given 
by 


ies sMM(A *tij—1, Aj), 


Uj—2 


u; = MM (ac 


where A*t; = At; ,1 — Ati;_1, and At; ,1 = %j41 —ti;. Here, the MinMod 
limiter (MM) is defined by 


min,{z,p} if x, > 0, 


MM(a1,22,...) = ¢ maxp{zp} if z, <0, 
0 otherwise. 


2.3 Essentially nonoscillatory reconstruction 


Spurious oscillations can appear in the numerical solution, if the big stencil 
which defines popt,;(a), contains a discontinuity or large gradients. To avoid 
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such a problem a WENO procedure is designed that smoothly adapts the 
stencil in the neighbourhood of the singularity. 


The principle of CWENO (central WENO) procedure which is defined 
in [8,24] is extended, and an ENO interpolant is constructed as a convex 
combination of polynomials that are based on different stencils; that is, 


R;(z) = S° wip} (2), w! > 0, yw =1, i€ {0,1,2,c}. (10) 


Here, R;(x) is the nonoscillatory reconstruction on J;. p(x) is a polynomial 
which is defined on stencil S, also the polynomials p}(x), pj (x), and p3(x) 
are computed in (f@). In order to simplify the notations, the upper index 
j will be omitted; remembering that the weights and the four polynomials 
change from cell to cell. Now, the polynomial p,(a) is computed such that 
the convex combination (0), will be fifth-order accurate in smooth regions. 
It must, therefore, satisfy the following equation: 


Popt (x) = Copo(x) +Cipi (x) + Cape2(x) +Cepe(x), yc = 17€ {0,1,2,c}, 


(11) 
where constants C; represent linear or ideal weights for (10). A straightfor- 
ward calculation shows that any symmetric choice of constants C; in (MH) 
provides the desired accuracy. This property must be contrasted with clas- 
sical upwind WENO schemes. Therefore, we make the choice Co = C2 = 
t Ci = i and C, = s. Then, the polynomial p,(a) is obtained explicitly as 
follows: 


Pe(£) = [Popt(x) — Copo(z) — Cipi(x) — Cope(x)]/Ce Va ed;. (12) 


In order to complete the reconstruction of R;(x), the ideal weights C;, 7 € 
{0, 1,2, c}, are changed to nonlinear weights w;, with objective of maintaining 
the fifth-order accuracy for smooth solutions and nonoscillatory performance 
near discontinuities. In smooth regions for achieving the optimal interpo- 
lation (1), the weights w; must smoothly converge to the ideal weights C; 
as h approaches zero. Also, in regions where a discontinuity does exist, the 
weights should effectively remove the contribution of stencils that contain the 
discontinuity. 


2.4 The nonoscillatory weights 


In order to fully determine the scheme, it is required to specify the nonoscil- 
latory weights,. Then, following [21,82], the nonlinear weights are written as 
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Ay C; ‘ 
i= , a= k € {0,1,2,c}. 13 
° a Ok _ (e+ 19;)?’ : { ch ve) 
The constants C;, 7 € {0,1,2,c}, in Eq. (3) ane chosen to be the same as 
in the equation (1); that is, Co = C2 = 3,C = j, and C, = 5. 


The smoothness indicators, IS;, are responsible for detecting large gradi- 
ents or discontinuities and automatically switch to the stencil that generates 
the least oscillatory reconstruction in such cases. The smoothness indicator 
is defined as [ZI] 


ks 


2 
> _ d ; 
IS; = h2*k 1 x | (ra) dx, LE {0, 1, 2,c}. 
k i 


The constant €, taken to prevent the denominator from vanishing, can range 
from 10~° to 107", [21]. As is said in [12] and [19], different values of € change 
the order of convergence of the scheme and yield different sharpness near 
discontinuities. In this paper, €e = 10~® is used. The smoothness indicators 
of smooth and discontinuous stencils yield 19, = O(h?) and IS, = O(1), 
respectively. A direct computation based on equations ({) and (J) yields: 


1 
IS, =—— (9tiy — 12054. 1) + 3tj4271-r) + (1 — 7) 96u4)? 


8100 
13 
3700 Oe — 66u;4(1- r) + 9U;420— —r) + (1 — r)48u;, 
781 7 : n ie 
162000 (72u,; 96u;4(1—r) + 24U542(1-r) + (1 = r)48u;) ; i! 0, 2, 
9 | 7 ca es 
IS ~ 9916 (tij-1 Uj+1 16u‘,)? + Ip (Hi-1 = 2u; + tij41)? 
449856 _ _ 
 Fa35q (Hit — Wit + 25)”, 
and 
#6 _ (23957 690709, 349134051 20591 ‘i 
¢ ~\37500 7+? ~ 126000 2+! * ga00 “2 126000 "7-1 126000 "9-274? 
1478263 __ 7423 68419 4051 _ ” 
126000 2+! 350 “7 * 210002"! — 126000 “9 -27"9+4 
1432397423, 34913 i; 
8400 2 350 2-1" 400 “97275 
1478263 _- 690709 Nien (23987 5 
126000 2-1 126000 “2~??"3-1 © 37500 9-27%5-2 
2464 , 649_ 2486 | 2486 _ 649 


U; Uj—2 + Uj-1 Uj+1 + Uij42)u; 
1125 7% 375 375 375 375 
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2.5 The new scheme for nonlinear degenerate parabolic 
PDEs 


From equations (@) and (9), it is required to approximate b(u), and vz. The 
final form of the left- and right-biased points are given by 


U4. = die {0,1,2,c} Wi X a) (7542), 


a » 4€{0,1,2,c} w, x ot) (x,;_1), 


Used + = Pole; 43) = (22d; + 9Uj41 — Tj42 + 8uj)/30, 


Ucn = Pale p44) = (Gj—1 + 15a; + Qj 41 + Buj)/18, 


Used | = Pol tj44) = (—2ij—2 + 13tij1 + 194, + 24u/)/30, 


Uist t= Pel%z41) 
= (30%; 2 — 2050;_1 + 291t, + 277%; 41 — 33%) — 176u',) /360, 
Uv. ? = po(x;—1) = (190, + 130541 = 20542 = 24v;,)/30, 
+1). _ (95. Sa... Rail 
Uy = pi(x;_1) = (20j;-1 + 150; + 0j41 — 8v;)/18, 


vt zi : = po(#j_1) = (—Oj-2 + 90j_1 + 220, — 8v4)/30, 


Wh) += palay-g) 
= (—330;_2 + 2770;-1 + 2910; — 2050,;+1 + 300j42 + 176v;)/360. 
Therefore, 
a Bt a) PW a) 
ha os (1a 
dt; (t j+5 j-%5 
ee) = ta 2 + S05, t,u,) = F(G;(t)), 


and the semidiscrete scheme (4) is discretized in time by a third order TVD 
RungeKutta method [B3], which is given by 
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3 Numerical results 


This section gives some numerical examples to confirm the accuracy, robust- 
ness, and high-resolution properties of the proposed method. Example 1 is 
three-dimensional heat equation with a known exact solution and is used 
for checking the accuracy and order of convergence [I5]. Example 2 is the 
two-dimensional porous medium equation [26]. In Example 3, a nonlinear 
problem is presented and suggested as test for multidimensional nonlinear 
convection-diffusion problems [8,22]. Next, in Example 4, a scalar two- 
dimensinal reaction-diffusion equation is considered [7]. Example 5 taken 
from [20,23], is presented to cover the case strongly degenerate parabolic 
equation. Finally, Example 6 is considered to solve a Turing model of bio- 
logical pattern formation [Z@. This example will demonstrate that the new 
technique also handles systems of degenerate equations. The constant CFL 
number c = 0.3 is considered. Then, the time step is determined by 


Also, in the following, CWENO5 is used to denote the new proposed method 
in the current paper. 


Example 1. For x = (a, y, z), S(x,t,u) = 0, and b;(u) =u, j = 1, 2,3, con- 
sider the three-dimensional linear initial-value problem with periodic bound- 
ary condition [5], 


Ut = Uga + Uyy F Uzz, X 2 = ( ra) be, 
u(x, y, z,0) = sin(x) sin(y) sin(z). 


The closed analytical form solution for this problem is 
u(x, y, z,t) = exp(—t#) sin(x) sin(y) sin(z). 


The error is computed in the discrete L,, and L, norms and, respectively, 
defined by 


I|ulloo = maxi, j,k [Uij,e1, 


lJul|a = Lad [ui jelAc?. 
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Table 1: Errors for Example 1 at T=2 


N Lo error L. order JL; error JL, order CPU 


10 0.192(-2) - 0.847(-3) = 08.97 
20 0.601(-4) 5.000 0.266(-4) 4.999 15.01 
40 0.189(-5) 4.999 —0.822(-6) 5.021 32.84 
80 0.579(-7) 5.030 ~—-0.257(-7) 5.000 ~—«6 8.10 
160 0.163(-8) 5.151 0.804(-9) 5.000 142.90 


Table 2: Errors and orders of convergence for Example 1 with different linear weights 
at T=2. Top:Co = Cy = Cz = Ce = +, Middle:Cyp = Cp = 4,C1 = 3,C2 = 3 
1 


16? 8? 
Bottom:Co + Ci =Co2 = Ce = 5- 


N Lo error L. order JL, error JL, order CPU 


10 (0.115(-2) - 0.788(-3) - 09.12 
20  0.400(-4) 4.845 0.295(-4) 4.739 15.15 
40  0.110(-5) 5.184 0.910(-6) 5.018 33.01 
80 0.347(-7) 4.986 0.290(-7) 4.971 70.11 
160  0.107(-8) 5.019 0.863(-9) 5.070 139.98 
10 -0.149(-2) - 0.102(-2) - 10.01 
20 = 0.453(-4) 5.039 0.341(-4) 4.902 14.93 
40  0.190(-5) 4.575 0.971(-6) 5.134 30.99 
80 —-0.338(-7) 5.812 0.281(-7) 5.110 73.12 
160  0.103(-8) 5.036 0.877(-9) 5.001 138.88 
10 -0.171(-2) - 0.145(-2) - 11.55 
20 = 0.148(-3) 3.530 0.139(-3) 3.382 16.01 
40  0.170(-4) 3.122 0.182(-4) 2.933 32.18 
80 = 0.125(-5) 3.765 0.201(-5) 3.178 69.87 
160  0.873(-7) 3.839 0.156(-6) 3.687 140.11 


In Table f, the respective errors and convergence rates by CWENO5 are 
given. The results verify the fifth-order accuracy both in the DL, and in the 
Ly, norms. The L; and LD. errors and orders of convergence by CWENO5 
with different linear weights are reported in Table 2 As is expected, any 
symmetric choice of constants C; enables the fifth-order accuracy. 


Example 2. Consider the porous medium equation (PME) as [26] 
uz, = A(u™), m> 1, (16) 


where u = u(x,t), x € Q C R*%. This equation describes various diffusion 
processes, such as the flow of an isentropic gas through a porous medium 
where u is the density of the gas required to be non-negative and u™~! is the 
pressure of the gas. The PME equation degenerates at points where u = 0, 
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Table 3: Errors and orders of convergence for Example 2 at T=2. Top to Bottom 
respectively: m= 2, m=4, m=6, andm=8. 


N Lo error L. order JL, error JL, order CPU 


50 =—-0.012(-1) - 0.429(-1) - 10.62 
100 0.010(-1) 0.263 0.172(-1) 1.318 19.88 
200 0.006(-1) 0.737 0.036(-1) 2.256 35.17 
400 0.002(-1) 1.585 0.006(-1) 2.585 78.41 
50 = -0.117(-1) - 0.248 - 13.51 
100 0.083(-1) 0.495 0.639(-1) 1.956 23.23 
200 0.069(-1) 0.266 0.256(-1) 1.319 37.11 
400  0.054(-1) 0.353 0.092(-1) 1.476 86.12 
50 =: 0.216 (-1) - 0.353 - 16.55 
100 0.232(-1) -0.103 0.135 1.386 27.51 
200 0.179(-1) 0.374 0.725(-1) 0.896 41.68 
400  0.134(-1) 0.417 0.248(-1) 1.547 95.90 
50 -0.337(-1) - 0.583 - 23.55 
100 -0.248(-1) 0.442 0.201 1.536 39.77 
200 0.213(-1) 0.219 0.752(-1) 1.418 53.13 
400 0.218(-1) -0.033 0.368(-1) 1.031 115.55 


resulting in the phenomenon of finite speed of propagation and sharp fronts. 
The classical solutions to the PME may not exist in general, even if the initial 
condition is smooth. The closed analytical form which is a weak solution for 
PME, is obtained by Zel’dovich and Kompaneets [88] and Barenblatt [5] 
in years 1950 and 1952, respectively. The solution of [5] has the following 
explicit form 


=. 
m—1 


u(x,t) =t-? (1 aixPe-*) | (17) 
+ 
where a = Seer k= lua and wu, = max(u,0). Now, we put d= 2; 


therefore, Barenblatt solution has no derivative at the points of the circle 
ga = 4 / Ce where a = +. The CWENO5 scheme is used to solve 


the PME (M3) for m = 2,4,6 and m = 8, where 2 = [—10,10]?. The initial 
condition is taken as the Barenblatt solution (I) at t = 1, and the boundary 
condition is chosen to be u = 0 on O02) for t > 1. 

Table B indicates the magnitude of errors obtained with CWENO5 scheme 
for different mesh sizes. Figure [shows the results of CWENO5 scheme at the 
Final time T = 3 with the computational domain 2 divided into 120 uniform 
cells. The CWENO5 scheme prevents the appearance of spurious solutions 
close to the circles. The contour plots of absolute difference between approx- 
imated solution and Barenblatt solution with m = 2,4 are shown in Figure [I] 
top right and bottom right, respectively. The CWENO5 scheme is sufficiently 
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Table 4: Errors and orders of convergence for Example 2 (m = 2) with different linear 
weights at T=2. Top:Co = C1 = Cz = Co = 4, Middle:Co = C2 = C1 31Ce 


16° 
Bottom:Co + Cy =C2 Co - 


8? 


N Lo error L,, order Ly, error JL, order CPU 


50 0.012(-1) ; 0.422(-1) ; 10.89 
100 0.009(-1) 0.415 0.165(-1) 1.354 20.85 
200 0.005(-1) 0.848  0.037(-1) 2.156 =~ 34.37 
400 0.002(-1) 1.322 ~=0.007(-1) ~—-2.402-~—80.51 
50 0.014(-1) : 0.418(-1) ‘ 11.11 
100 0.008(-1) 0.807 0.167(-1) 1.324 22.20 
200 0.006(-1) 0.415  —0.039(-1) 2.098 ~—-33.01 
400 0.003(-1) 1.000 ~=—-0.009(-1) «2.116 ~—81.52 
50 0.222(-1) : 0.338 : 09.95 
100 0.218(-1) 0.026 0.165 1.035 22.51 
200 0.157(-1) 0.473 ~—-0.722(-1) —s«:1.192-—-35.68 
400 0.131(-1) 0.261 ~=—-0.243(-1) «1.571 79.95 


0.6 


10. —5 r) 5 10 
10 0.06 
- 0.05 
0.04 
2 0.03 
5 0.02 
0.01 
-10 
-10 -5 oO 5 10 


Figure 1: Example 2 (two-dimensional porous medium equation). Top left: PME with 
m =2and N = 120, Top right: 30 contour lines of |wBarenblatt — U| with m = 2, Bottom 
left: PME with m = 4 and N = 120, Bottom right: 30 contour lines of |uBarenblatt — U| 
with m = 4. 
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Figure 2: Example 3 (two-dimensional Buckley—Leverett-type equation). Left: numerical 
solution at T= 0.5 with N = 40, Right: numerical solution at T = 1.0 with N = 40. 


accurate in the smooth regions, while on the nonsmooth regions absolute dif- 
ference with analytical form solution of Barenblatt is of order 107° and 107? 
for m = 2 and m = 4, respectively. 


Example 3. This example solves the Buckley—Leverett-type problem [B,22] 


with « = 0.1, the flux function of the form 

2 
u2 + (1— 1)?’ 
f(u) = g(u)(1 —5(1 — u)?). 


The initial function for this example is 


g(u) = 


1, (x — 0.25)? + (y — 0.25)? <5, 
ule) = i otherwise. 
This example includes gravitational effects in the x-direction. The numerical 
results at T = 0.5 and 1.0 are shown in Figure 2, with the computational 
domain [—3, 3]? divided into 40 x 40 uniform cells. The results compare well 
with those reported in [B]. Table Hl gives the respective errors and convergence 
rates by CWENOS. The numerical solutions compared with a “reference so- 
lution”, obtained by WENO6 [26] with N = 400. 


Example 4. This example solves the following initial-boundary-value prob- 
lem for a scalar reaction-diffusion equation [f/], where x = (x,y), 


uz = Ab(u) + S(x,t, u), (18) 


u(x, 0) = 0.5(1 + sin(1.1(a — cos(0.7y)))) cos(0.5(y — sin(1.32))), 
Vb(u)-n=0 on OQ x (0,7). (19) 
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Table 5: Errors and orders of convergence for Example 3 at T=0.5 


N Lo error L. order JL; error JL, order CPU 


40 0.228(-3) : 0.716(-4) : 33.63 
80 0.680(-5) 4.987  0.248(-5) 4.852 58.12 
160 0.215(-6) 4.998 0.799(-7) 4.955 ~—- 109.87 
320 0.671(-8) 4.992 0.239(-8) 5.060 ~—-219.23 


This problem may serve as a scalar prototype degenerate reaction-diffusion 
model. Here, the zero-flux boundary condition (1) implies that the reaction- 
diffusion domain is isolated from the external environment. For S(x, t,u) = 
Sw), (28) appears in [28] in an ecological setting, where u denotes the pop- 
ulation density of a species, and S(u) is its dynamics, where it is assumed 
that S(0) = 0 and $’(0) 4 0. For example, S(u) = u(1 — u) — u?/(1 + u?) 
corresponds to the population dynamics of the spruce band-worm [28] and 
models the growth of the population by a logistic expression and the rate of 
mortality due to predation by other animals. Authors of [Z| modified this 
expression by a radial spatial factor, and used 


S(x,t,u) = 10(exp(—5r)u(1 — u) + (exp(—5r) la): 


r= V(x — 0.5)? + (y — 0.5)?, 


which means that the birth of individuals is concentrated near the centre 
(0.5,0.5), and mortality increases with increasing distance from the centre. 
Most standard spatial models of population dynamics simply assume that 
b(u) = Du, where the constant diffusion coefficient D > 0 measures the 
dispersal efficiency of the species under consideration. In this example, the 
strongly degenerate diffusion coefficient [4,7] is used 


b(u) _ 0, uU< Ue, 
~ | Diu ue), otherwise, 


where u, > 0 is an assumed critical value of u beyond which diffusion will take 
place. Here, we set D = 1 and ue = $. The difficulty in the well-posedness 
analysis of the problem (TX) lies in the boundary condition (19) when b is 
strongly degenerate. It is quite difficult to give a correct formulation of the 
zero-flux boundary conditions. The computational domain [0,1]? is divided 
into 80 x 80 uniform cells. Figure 8 shows the numerical approximation at 
different time levels. The results compare well with those reported in [i]. 


Example 5. In this example, we focus on strongly degenerate parabolic 
problem. Consider the Burgers-like equation 


ue t+ (Ua + (U?)y = €(v(u)ur)e + €(Y(U)Uy)y. 


WENO schemes for multidimensional nonlinear degenerate... 


1 
0.6 
0.5 
0.4 
0.3 
0.2 
0.1 
ro) 
) 1 
1 0.5 
0.4 
0.3 
0.2 
0.1 
) 
) 1 


Figure 3: Example 4 (single-species reaction-diffusion model). From top left to bottom 
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right, it shows the numerical solution at times T’ = 0.0, 0.5, 1.0, 3.0. 


Table 6: Errors and orders of convergence for example 5 at T=0.5 


N Lo error L. order JL, error Ly, order CPU 
40  0.223(-4) - 0.118(-4) - 30.73 
80  0.715(-6) 4.963 0.362(-6) 5.000 53.98 
160 0.214(-7) 5.059 0.112(-7) 5.014 100.11 
320 0.599(-9) 5.162 0.334(-9) 5.054 213.23 
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0.55 


Figure 4: Example 5 (two-dimensional strongly degenerate parabolic problem). Top: 
numerical solution at T = 0.5 with N = 80, Bottom: one-dimensional cut along the y 
axis with different mesh sizes. 


Here, € = 0.1 is considered and 


(u) = £° |u| < 0.25, 
w=) 1 lal > 0.25. 


In this example, zero boundary conditions are applied. The initial condition 
is equal to —1 and 1 inside two circles of radius 0.4 centred at (0.5,0.5) and 
(—0.5, —0.5), respectively, and is zero elsewhere inside the square [—1.5, 1.5]?. 
The LD. and Ly errors of CWENO5 scheme are given in Table @ for different 
mesh sizes. The computational domain is divided into 80 x 80 uniform cells. 
The numerical results obtained by CWENO5 scheme are presented in Figure 
A, which compare well with those reported in [23,26]. The “reference solution” 
was obtained by WENO6 [26] scheme with N = 400. 


Example 6. Turing [84] wrote his paper “The chemical basis of morpho- 
genesis” in which he suggested an new theory. He hypothesized that the 
patterns, observed during embryonic development, occur in reaction to a spa- 
tial prepattern in biochemicals, which he termed morphogens. Cells would 
then answer to this prepattern by differentiating in a threshold-dependent 
way. Thus, Turing hypothesized that the patterns, observed in nature, such 
as pigmentation in animals, branching in trees, and skeletal structures are 
reflections of inhomogeneities in underlying biochemical signalling. 


The system he considered took the form 


ou > 
an DV~u+ f(u). 
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Figure 5: Example 6 (Turing’s model for biological pattern formation). Left: wu after 
1113 simulated time units, Right: v after 1113 simulated time units. 


Dy 

0 D, 
constant diffusion coefficients, and the nonlinear function f(u) is the reaction 
kinetics. Turing analysed the special case of equation (8) when there are 
two interacting morphogens, u = (u,v)’. They spread at rates D,, and 
D,, respectively, and they also undergo reaction kinetics defined by f(u) = 
(f(u,v),g(u,v))?. There are different forms of the reaction kinetics such as 
Gierer—Meinhardt [14] model 


Here, u is a vector of chemical condensations, D = | is a matrix 


u2 


os 2 
(+ ku)v’ g(u, v) = C4U — C5U 


f(u,v) =a — cqQu+cs3 


and the Schnakenberg [BT] model 


f(u,v) =e, —c_yute3u?v, — g(u,v) = cz — cu. 

The parameters {c_1,¢1,C2,c3} represent the deterministic rates of the re- 
actions. Now, consider the Schnakenberg model on the unit square domain 
Q = [0,1] x [0,1] with zero-flux boundary conditions. Let the parameter 
values be determined by D, = 1, Dy = 40, cy = 0.1, co = 0.9,c_1 = 1, 
and cz = 1. The initial conditions are random perturbations about the ho- 
mogeneous steady state. This system has a uniform positive steady state 
(u°, v°), where uw? = cy +c and v° = c2/(c; + c2)?. Figure H shows the for- 
mation of spot patterns. In this figure, the concentration values of u and v 
are illustrated which compare well with those reported in [27]. 


4 Conclusion 


When applying the standard WENO idea, as Liu, Shu, and Zhang did in [26], 
in some cases (odd numbers of nodes in small stencils) the linear weights 
don’t exist, and when they do exist, at least one of them is negative. In 
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this research, another strategy was used to circumvent these problems. This 
strategy consists in applying the idea of the local discontinuous Galerkin 
method. First, the degenerate parabolic equation was considered as a non- 
linear system of first order equations, and then this system was solved us- 
ing a fifth-order finite difference weighted essentially nonoscillatory (WENO) 
method for conservation laws. The Symmetrical WENO procedure of [I] for 
solving this system was used. It was observed that this procedure generates a 
fifth-order method in smooth regions and remains essentially nonoscillatory 
near discontinuities. Numerical examples showed that the new scheme can 
obtain fifth-order accuracy in smooth regions and prevent the appearance of 
spurious solutions close to discontinuities. 
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